﻿#pragma once

#include <iostream>
#include <optional>
#include <random>

#include <xy/math/math.h>

namespace xy
{

namespace geom
{

template <std::size_t M, std::size_t N, typename T> struct mat
{
  protected:
    T m_data[M * N];

  public:
    /**
     * @brief 基本空构造
     *
     */
    mat(T t = T{})
    {
        for (std::size_t i = 0; i < M * N; i++)
            m_data[i] = t;
    }

    /**
     * @brief 变长构造
     *
     * @tparam Ts
     * @return requires(... == N && (std::is_same_v<Ts, T> && ...)) vec(const Ts&... ts)
     */
    template <class... Ts>
        requires(sizeof...(Ts) == M * N && (std::is_same_v<Ts, T> && ...))
    mat(const Ts &...ts) : m_data{ts...}
    {
    }

    /**
     * @brief 拷贝构造
     *
     * @param[in] rhs
     */
    mat(const mat &rhs)
    {
        std::copy(rhs.data(), rhs.data() + M * N, m_data);
    }

    /**
     * @brief 拷贝赋值
     *
     * @param[in] rhss
     * @return mat&
     */
    const mat &operator=(const mat &rhs)
    {
        if (this != &rhs) [[likely]]
            std::copy(rhs.data(), rhs.data() + M * N, m_data);
        return *this;
    }

    /**
     * @brief 取反
     *
     * @return mat
     */
    mat operator-() const
    {
        mat ret;
        for (std::size_t i = 0; i < M * N; i++)
            ret.data()[i] = -m_data[i];
        return ret;
    }

    /**
     * @brief 获得元素
     *
     * @param[in] i
     * @param[in] j
     * @return T
     */
    T operator()(std::size_t i, std::size_t j) const
    {
        if (i >= M || j >= N)
            throw std::out_of_range("Index out of range");
        return m_data[i * N + j];
    }

    /**
     * @brief 获得元素
     *
     * @param[in] i
     * @param[in] j
     * @return T&
     */
    T &operator()(std::size_t i, std::size_t j)
    {
        if (i >= M || j >= N)
            throw std::out_of_range("Index out of range");
        return m_data[i * N + j];
    }

#define MAT_BINARY_OP(func, op)                                                                                        \
    mat func(T rhs) const                                                                                              \
    {                                                                                                                  \
        mat ret;                                                                                                       \
        for (std::size_t i = 0; i < M; i++)                                                                            \
            for (std::size_t j = 0; j < N; j++)                                                                        \
                ret(i, j) = m_data[i * N + j] op rhs;                                                                  \
        return ret;                                                                                                    \
    }

    MAT_BINARY_OP(operator*, *)
    MAT_BINARY_OP(operator/, /)
#undef MAT_BINARY_OP

#define MAT_BINARY_OP(func, op)                                                                                        \
    mat func(const mat &rhs) const                                                                                     \
    {                                                                                                                  \
        mat ret;                                                                                                       \
        for (std::size_t i = 0; i < M; i++)                                                                            \
            for (std::size_t j = 0; j < N; j++)                                                                        \
                ret(i, j) = m_data[i * N + j] op rhs(i, j);                                                            \
        return ret;                                                                                                    \
    }

    MAT_BINARY_OP(operator+, +)
    MAT_BINARY_OP(operator-, -)
#undef MAT_BINARY_OP

#define MAT_BINARY_OP(func, op)                                                                                        \
    mat &func(T rhs)                                                                                                   \
    {                                                                                                                  \
        for (std::size_t i = 0; i < M; i++)                                                                            \
            for (std::size_t j = 0; j < N; j++)                                                                        \
                m_data[i * N + j] op rhs;                                                                              \
        return *this;                                                                                                  \
    }

    MAT_BINARY_OP(operator*=, *=)
    MAT_BINARY_OP(operator/=, /=)
#undef MAT_BINARY_OP

#define MAT_BINARY_OP(func, op)                                                                                        \
    const mat &func(const mat &rhs)                                                                                    \
    {                                                                                                                  \
        for (std::size_t i = 0; i < M; i++)                                                                            \
            for (std::size_t j = 0; j < N; j++)                                                                        \
                m_data[i * N + j] op rhs(i, j);                                                                        \
        return *this;                                                                                                  \
    }

    MAT_BINARY_OP(operator+=, +=)
    MAT_BINARY_OP(operator-=, -=)
#undef MAT_BINARY_OP

    /**
     * @brief 矩阵乘法
     *
     * @tparam P
     * @param[in] mat
     * @return mat<M, P, T>
     */
    template <std::size_t P> mat<M, P, T> operator*(const mat<N, P, T> &matrix) const
    {
        T ret[M * P], data[N * P];
        std::fill(ret, ret + M * P, T{});
        std::copy(matrix.data(), matrix.data() + N * P, data);
        for (std::size_t i = 0; i < M; i++)
            for (std::size_t j = 0; j < P; j++)
                for (std::size_t k = 0; k < N; k++)
                    ret[i * P + j] += m_data[i * N + k] * data[k * P + j];

        mat<M, P, T> result;
        std::copy(ret, ret + M * P, result.data());
        return result;
    }

    /**
     * @brief 指标返回
     *
     * @param[in] i
     * @return T
     */
    T operator[](std::size_t i) const
    {
        if (i >= M * N)
            throw std::out_of_range("Index out of range");
        return m_data[i];
    }

    /**
     * @brief 指标返回
     *
     * @param[in] i
     * @return T&
     */
    T &operator[](std::size_t i)
    {
        if (i >= M * N)
            throw std::out_of_range("Index out of range");
        return m_data[i];
    }

    /**
     * @brief 转置
     *
     * @return mat<N, M, T>
     */
    mat<N, M, T> transposed() const
    {
        mat<N, M, T> ret;
        for (std::size_t i = 0; i < M; i++)
            for (std::size_t j = 0; j < N; j++)
                ret(j, i) = m_data[i * N + j];
        return ret;
    }

    /**
     * @brief 取模
     *
     * @return T
     */
    T modulus() const
    {
        return std::sqrt(square_modulus());
    }

    /**
     * @brief 模平方
     *
     * @return T
     */
    T square_modulus() const
    {
        T r{};
        for (std::size_t i = 0; i < M * N; i++)
            r += m_data[i] * m_data[i];
        return r;
    }

    /**
     * @brief 正则化.
     * @see
     */
    mat &normalize()
    {
        T mod = modulus();
        if (mod < math::numeric::distance)
            throw std::runtime_error("Divide by zero");
        for (std::size_t i = 0; i < M * N; i++)
            m_data[i] /= mod;
        return *this;
    }

    /**
     * @brief 返回正则化
     *
     * @return mat
     */
    mat normalized() const
    {
        T mod = modulus();
        if (mod < math::numeric::distance)
            throw std::runtime_error("Divide by zero");
        mat ret;
        for (std::size_t i = 0; i < M * N; i++)
            ret.data()[i] = m_data[i] / mod;
        return ret;
    }

    /**
     * @brief 绝对值
     *
     * @return mat&
     */
    mat &absolute()
    {
        for (std::size_t i = 0; i < M * N; i++)
            m_data[i] = std::abs(m_data[i]);
        return *this;
    }

    /**
     * @brief 返回绝对值向量
     *
     * @return vec
     */
    mat absoluted() const
    {
        mat ret;
        for (std::size_t i = 0; i < M * N; i++)
            ret[i] = std::abs(m_data[i]);
        return ret;
    }

    /**
     * @brief 到点的距离
     *
     * @param[in] mat
     * @return T
     */
    T distance(const mat &mat) const
    {
        return std::sqrt(square_distance(mat));
    }

    /**
     * @brief 到点的距离平方
     *
     * @param[in] mat
     * @return T
     */
    T square_distance(const mat &mat) const
    {
        return (*this - mat).square_modulus();
    }

    /**
     * @brief 比较运算
     *
     * @param[in] mat
     * @return true
     * @return false
     */
    bool operator==(const mat &mat) const
    {
        return distance(mat) < math::numeric::half;
    }

#define MAT_VAL_OP(func, i)                                                                                            \
    T func() const                                                                                                     \
    {                                                                                                                  \
        static_assert(i < M * N, "Index out of range");                                                                \
        return m_data[i];                                                                                              \
    }                                                                                                                  \
    T &func()                                                                                                          \
    {                                                                                                                  \
        static_assert(i < M * N, "Index out of range");                                                                \
        return m_data[i];                                                                                              \
    }

    MAT_VAL_OP(x, 0)
    MAT_VAL_OP(y, 1)
    MAT_VAL_OP(z, 2)
    MAT_VAL_OP(w, 3)
    MAT_VAL_OP(r, 0)
    MAT_VAL_OP(g, 1)
    MAT_VAL_OP(b, 2)
    MAT_VAL_OP(a, 3)
    MAT_VAL_OP(front, 0)
    MAT_VAL_OP(back, M *N - 1)
#undef MAT_VAL_OP

    /**
     * @brief 获得数据数组
     *
     * @return const T*
     */
    const T *data() const noexcept
    {
        return m_data;
    }

    /**
     * @brief 获得数据数组
     *
     * @return T*
     */
    T *data() noexcept
    {
        return m_data;
    }

    /**
     * @brief 单位阵
     *
     * @return mat
     */
    static mat identity()
    {
        static_assert(M == N, "Identity matx must be square");
        mat ret;
        for (std::size_t i = 0; i < M; ++i)
            ret(i, i) = 1;
        return ret;
    }
};

// 高精度版本
using mat4x4 = mat<4, 4, double>;
using mat3x3 = mat<3, 3, double>;
using mat2x2 = mat<2, 2, double>;

// 浮点版本
using mat4x4f = mat<4, 4, float>;
using mat3x3f = mat<3, 3, float>;
using mat2x2f = mat<2, 2, float>;

template <std::size_t N, typename T> using vec = mat<N, 1, T>;

// 高精度版本
using vec2 = vec<2, double>;
using vec3 = vec<3, double>;
using vec4 = vec<4, double>;

// 浮点版本
using vec2f = vec<2, float>;
using vec3f = vec<3, float>;
using vec4f = vec<4, float>;

// 整型版本
using vec2i = vec<2, int>;
using vec3i = vec<3, int>;
using vec4i = vec<4, int>;

/**
 * @brief 输出函数
 *
 * @tparam M
 * @tparam N
 * @tparam T
 * @param os
 * @param m
 * @return std::ostream&
 */
template <std::size_t M, std::size_t N, class T>
inline std::ostream &operator<<(std::ostream &os, const mat<M, N, T> &m)
{
    static constexpr std::size_t row = N == 1 ? 1 : M;
    static constexpr std::size_t col = N == 1 ? M : N;

    os << "{ ";
    if constexpr (row > 1)
        std::cout << std::endl;

    for (std::size_t i = 0; i < row; i++)
    {
        for (std::size_t j = 0; j < col; j++)
            os << m[i * col + j] << " ";

        if constexpr (row > 1)
            std::cout << std::endl;
    }
    os << "}";
    return os;
}

// 定义向量运算
#define MAT_OP(func, op)                                                                                               \
    template <std::size_t N, class T> inline vec<N, T> func(const vec<N, T> &v1, const vec<N, T> &v2)                  \
    {                                                                                                                  \
        vec<N, T> ret;                                                                                                 \
        for (std::size_t i = 0; i < N; i++)                                                                            \
            ret[i] = v1[i] op v2[i];                                                                                   \
        return ret;                                                                                                    \
    }

MAT_OP(operator*, *)
MAT_OP(operator/, /)

#undef MAT_OP

#define MAT_OP(func, op)                                                                                               \
    template <std::size_t N, class T> inline vec<N, T> &func(vec<N, T> &v1, const vec<N, T> &v2)                       \
    {                                                                                                                  \
        for (std::size_t i = 0; i < N; i++)                                                                            \
            v1[i] op v2[i];                                                                                            \
        return v1;                                                                                                     \
    }

MAT_OP(operator*=, *=)
MAT_OP(operator/=, /=)

#undef MAT_OP

/**
 * @brief 类型转换
 *
 * @tparam N
 * @tparam U
 * @tparam V
 * @param[in] v
 * @return vec<N, U>
 */
template <class U, std::size_t N, class V> inline vec<N, U> convert(const vec<N, V> &v)
{
    vec<N, U> ret;
    for (std::size_t i = 0; i < N; i++)
        ret[i] = static_cast<U>(v[i]);
    return ret;
}

/**
 * @brief 内积运算
 *
 * @tparam N
 * @tparam T
 * @param[in] v1
 * @param[in] v2
 * @return T
 */
template <std::size_t N, class T> inline T dot(const vec<N, T> &v1, const vec<N, T> &v2)
{
    T r{};
    for (std::size_t i = 0; i < N; i++)
        r += v1[i] * v2[i];
    return r;
}

/**
 * @brief 叉积
 *
 * @tparam T
 * @param[in] v1
 * @param[in] v2
 * @return vec<3, T>
 */
template <class T> inline vec<3, T> cross(const vec<3, T> &v1, const vec<3, T> &v2)
{
    return {v1.y() * v2.z() - v1.z() * v2.y(), v1.z() * v2.x() - v1.x() * v2.z(), v1.x() * v2.y() - v1.y() * v2.x()};
}

/**
 * @brief 计算夹角弧度
 *
 * @tparam N
 * @tparam T
 * @param v1
 * @param v2
 * @return T
 */
template <std::size_t N, class T> inline T angle(const vec<N, T> &v1, const vec<N, T> &v2)
{
    return std::acos(math::clamp(dot(v1, v2) / (v1.modulus() * v2.modulus()), T{-1}, T{1}));
}

/**
 * @brief 返回正交向量
 *
 * @tparam N
 * @tparam T
 * @param[in] v
 * @return vec<N, T>
 */
template <std::size_t N, class T> inline vec<N, T> orthogonal_vec(const vec<N, T> &v)
{
    // 巨坑：有些版本编译器不支持将 static_assert(false) 写在 else 里
    static_assert(N <= 3, "orthogonal_vec is not defined for N > 3");
    if constexpr (N == 2)
        return {-v.y(), v.x()};
    else if constexpr (N == 3)
    {
        if (std::abs(v.x()) < math::numeric::distance)
            return {T{}, -v.z(), v.y()};
        else if (std::abs(v.y()) < math::numeric::distance)
            return {v.z(), T{}, -v.x()};
        else if (std::abs(v.z()) < math::numeric::distance)
            return {-v.y(), v.x(), T{}};
        return {T{}, -v.z(), v.y()};
    }
}

/**
 * @brief 透视矫正，注意需要 screen_mapping 将 w 取倒数
 *
 * @tparam T
 * @param[in] coords
 * @param[in] z0
 * @param[in] z1
 * @param[in] z2
 * @return vec<3, T>
 */
template <class T> inline vec<3, T> perspective_correct(const vec<3, T> &coords, T w0, T w1, T w2)
{
    vec<3, T> ret = coords;
    ret[0] *= w0;
    ret[1] *= w1;
    ret[2] *= w2;
    T w = T{1} / (ret[0] + ret[1] + ret[2]);
    ret[0] *= w;
    ret[1] *= w;
    ret[2] *= w;
    return ret;
}

/**
 * @brief 笛卡尔坐标转齐次坐标
 *
 * @tparam N
 * @tparam T
 * @param[in] v
 * @param[in] w
 * @return vec<N + 1, T>
 */
template <std::size_t N, class T> inline vec<N + 1, T> cartesian_to_homogeneous(const vec<N, T> &v, T w = 1)
{
    vec<N + 1, T> p;
    std::copy(v.data(), v.data() + N, p.data());
    p[N] = w;
    return p;
}

/**
 * @brief 齐次坐标转笛卡尔坐标
 *
 * @tparam N
 * @tparam T
 * @param[in] p
 * @return vec<N, T>
 */
template <std::size_t N, class T> inline vec<N - 1, T> homogeneous_to_cartesian(const vec<N, T> &p)
{
    static_assert(N > 1, "homogeneous_to_cartesian is not defined for N < 2");
    vec<N - 1, T> v;
    std::copy(p.data(), p.data() + N - 1, v.data());
    return v / p[N - 1];
}

/**
 * @brief 生成 [x0,x1] x [x0,x1] x ... 方形区域内的随机点
 *
 * @tparam N
 * @tparam T
 * @param[in] x0
 * @param[in] x1
 * @return vec<N, T>
 */
template <std::size_t N, class T> inline vec<N, T> random_point(T x0, T x1)
{
    static std::mt19937 generator(std::random_device{}());
    if constexpr (std::is_floating_point_v<T>)
    {
        static std::uniform_real_distribution<T> uniform(x0, x1);
        vec<N, T> p;
        std::generate(p.data(), p.data() + N, [&]() { return uniform(generator); });
        return p;
    }
    else
    {
        static std::uniform_int_distribution<T> uniform(x0, x1);
        vec<N, T> p;
        std::generate(p.data(), p.data() + N, [&]() { return uniform(generator); });
        return p;
    }
}

/**
 * @brief 求逆
 *
 * @return std::optional<mat<N, N, T>>
 */
template <std::size_t N, class T> inline std::optional<mat<N, N, T>> invert(const mat<N, N, T> &matrix)
{
    // 初始化单位阵
    auto m = mat<N, N, T>::identity();
    auto old = matrix;

    T max;             // 最大值
    std::size_t index; // 最大值位置

    for (std::size_t i = 0; i < N; i++)
    {
        max = T{};
        index = i;

        // 寻找最大列元作为主元
        for (std::size_t j = i; j < N; j++)
        {
            T temp = std::abs(old(j, i));
            if (temp > max)
            {
                max = temp;
                index = j;
            }
        }

        if (max < math::numeric::distance)
            return std::nullopt;

        // 交换主元
        for (std::size_t q = 0; q < N; q++)
        {
            std::swap(old(index, q), old(i, q));
            std::swap(m(index, q), m(i, q));
        }

        // Gauss 消去
        for (std::size_t j = 0; j < N; j++)
        {
            T temp = old(i, i);
            if (i == j)
            {
                for (std::size_t k = 0; k < N; k++)
                {
                    m(j, k) /= temp;
                    old(j, k) /= temp;
                }
                continue;
            }
            temp = old(j, i) / old(i, i);
            for (std::size_t k = 0; k < N; k++)
            {
                m(j, k) -= temp * m(i, k);
                old(j, k) -= temp * old(i, k);
            }
        }
    }
    return m;
}

/**
 * @brief 平移矩阵
 *
 * @tparam T
 * @param[in] d
 * @return mat<4, 4, T>
 */
template <typename T> inline mat<4, 4, T> translate(const mat<3, 1, T> &d)
{
    auto A = mat<4, 4, T>::identity();
    A(0, 3) = d.x();
    A(1, 3) = d.y();
    A(2, 3) = d.z();
    return A;
}

/**
 * @brief 伸缩矩阵
 *
 * @tparam T
 * @param[in] scale
 * @return mat<4, 4, T>
 */
template <typename T> inline mat<4, 4, T> scale(const mat<3, 1, T> &scale)
{
    auto S = mat<4, 4, T>::identity();
    S(0, 0) = scale.x();
    S(1, 1) = scale.y();
    S(2, 2) = scale.z();
    return S;
}

/**
 * @brief 旋转矩阵
 *
 * @tparam T
 * @param[in] axis
 * @param[in] angle
 * @return mat<4, 4, T>
 */
template <typename T> inline mat<4, 4, T> rotate(const mat<3, 1, T> &axis, T angle)
{
    auto kx = axis.x(), ky = axis.y(), kz = axis.z();
    auto c = std::cos(angle), s = std::sin(angle);

    mat<4, 4, T> R;
    R(0, 0) = kx * kx * (1 - c) + c;
    R(0, 1) = kx * ky * (1 - c) - kz * s;
    R(0, 2) = kx * kz * (1 - c) + ky * s;
    R(1, 0) = kx * ky * (1 - c) + kz * s;
    R(1, 1) = ky * ky * (1 - c) + c;
    R(1, 2) = ky * kz * (1 - c) - kx * s;
    R(2, 0) = kx * kz * (1 - c) - ky * s;
    R(2, 1) = ky * kz * (1 - c) + kx * s;
    R(2, 2) = kz * kz * (1 - c) + c;
    R(3, 3) = 1;
    return R;
}

/**
 * @brief 计算正交投影矩阵，归一化映射到 [-1,1]^3
 * 输入参数 zNear, zFar 取正，因此 z 相关参数取负
 *
 * @tparam T
 * @param[in] left
 * @param[in] right
 * @param[in] bottom
 * @param[in] top
 * @param[in] zNear
 * @param[in] zFar
 * @return mat<4, 4, T>
 */
template <typename T> inline mat<4, 4, T> ortho(T left, T right, T bottom, T top, T zNear, T zFar)
{
    auto A = mat<4, 4, T>::identity();

    A(0, 0) = 2 / (right - left);
    A(1, 1) = 2 / (top - bottom);
    A(2, 2) = -2 / (zFar - zNear);

    A(0, 3) = -(left + right) / (right - left);
    A(1, 3) = -(top + bottom) / (top - bottom);
    A(2, 3) = -(zFar + zNear) / (zFar - zNear);
    // z 分量的分子分母都取负，因此抵消

    return A;
}

/**
 * @brief 计算透视投影矩阵
 * 透视矩阵
 * n: 近平面
 * f: 远平面
 * n, f 与 z 同号
 * 由于相机朝向 -z，输入参数 zNear, zFar 取正，因此矩阵需要取负
 * A = ( -n,  0,        0,  0 )
 *     (  0, -n,        0,  0 )
 *     (  0,  0, -(n + f), nf )
 *     (  0,  0,       -1,  0 )
 * n = -zNear, f = -zFar
 * A = ( zNear,     0,             0,            0 )
 *     (     0, zNear,             0,            0 )
 *     (     0,     0,  zNear + zFar, zNear * zFar )
 *     (     0,     0,            -1,            0 )
 *
 * @tparam T
 * @param[in] fovy 视角度数
 * @param[in] aspect 横纵比
 * @param[in] zNear 最近面
 * @param[in] zFar 最远面
 * @return mat<4, 4, T>
 */
template <typename T> inline mat<4, 4, T> perspective(T fovy, T aspect, T zNear, T zFar)
{
    auto A = mat<4, 4, T>::identity();

    // 将平截头体映射为长方体
    A(0, 0) = zNear;
    A(1, 1) = zNear;
    A(2, 2) = zNear + zFar;
    A(2, 3) = zNear * zFar;
    A(3, 2) = -1;
    A(3, 3) = 0;

    // 做正交投影
    T bottom = -zNear * std::tan(math::degrees_to_radians(fovy / 2));
    T top = -bottom;
    T left = bottom * aspect;
    T right = -left;

    auto O = ortho(left, right, bottom, top, zNear, zFar);
    return O * A;
}

/**
 * @brief 计算视图矩阵
 * 透视变换
 * g: 视线方向
 * t: 视线上方
 * (g x t, t, g) -> (e1, e2, -e3)
 * 逆变换 (e1, e2, e3) -> (g x t, t, -g)
 * R^{-1} = ( x_(g x t), x_t, x_(-g), 0 )
 *          ( y_(g x t), y_t, y_(-g), 0 )
 *          ( z_(g x t), z_t, z_(-g), 0 )
 *          (         0,   0,      0, 1 )
 * R^{-1}e1 = g x t, R^{-1}e2 = t, R^{-1}e3 = -g
 * R = (R^{-1})^T
 *
 * @tparam T
 * @param[in] pos 视点
 * @param[in] origin 视点指向
 * @param[in] up 视线方向
 * @return mat<4, 4, T>
 */
template <typename T>
inline mat<4, 4, T> look_at(const mat<3, 1, T> &pos, const mat<3, 1, T> &origin, const mat<3, 1, T> &up)
{
    auto Tr = mat<4, 4, T>::identity();
    Tr(0, 3) = -pos.x();
    Tr(1, 3) = -pos.y();
    Tr(2, 3) = -pos.z();

    auto g = (origin - pos).normalize();
    auto x = cross(g, up).normalize();
    auto t = cross(x, g);

    auto R = mat<4, 4, T>::identity();
    R(0, 0) = x.x();
    R(0, 1) = x.y();
    R(0, 2) = x.z();

    R(1, 0) = t.x();
    R(1, 1) = t.y();
    R(1, 2) = t.z();

    R(2, 0) = -g.x();
    R(2, 1) = -g.y();
    R(2, 2) = -g.z();

    return R * Tr;
}

/**
 * @brief 屏幕坐标映射
 *
 * @param[in] v
 * @param[in] width
 * @param[in] height
 * @return geom::vec4f
 */
template <typename T> inline vec<4, T> screen_mapping(const vec<4, T> &v, int width, int height)
{
    T w = v.w();
    auto p = v / w;
    T x = (p.x() + T{1}) * (width - 1) / T{2};
    T y = -(p.y() - T{1}) * (height - 1) / T{2};
    T z = p.z();

    // 透视矫正 1 / w
    return {x, y, z, T{1} / w};
}

/**
 * @brief 深度映射，解决深度分布问题
 * 由于透视矩阵除以 z 坐标，因此 z 值与 1/z 成正比
 *
 * @tparam T
 * @param z
 * @param zNear
 * @param zFar
 * @return std::uint32_t
 */
template <typename T> inline std::uint32_t depth_mapping(T z, T zNear, T zFar)
{
    // 注意输入 z 是负数，而 zNear, zFar 取正
    z = T{1} / (z * (T{1} / zFar - T{1} / zNear) + T{1} / zNear);
    z = (z - zNear) / (zFar - zNear);
    // 0.9 防止溢出，使用 1.0 乘法会得到 0
    return 0xffffffff * math::clamp(z, T{0}, T{0.9});
}

/**
 * @brief 计算反射向量
 *
 * @tparam N
 * @tparam T
 * @param[in] v
 * @param[in] n
 * @return vec<N, T>
 */
template <int N, class T> inline vec<N, T> reflect(const vec<N, T> &v, const vec<N, T> &n)
{
    return v - n * dot(v, n) * T{2};
}

/**
 * @brief 计算折射向量
 *
 * @tparam N
 * @tparam T
 * @param[in] v
 * @param[in] n
 * @param[in] eta = eta / eta'
 * @return std::optional<vec<N, T>>
 */
template <int N, class T> inline std::optional<vec<N, T>> refract(const vec<N, T> &v, const vec<N, T> &n, T eta)
{
    auto vn = dot(v, n);
    auto delta = T{1} - (eta * eta) * (T{1} - (vn * vn));
    if (delta < 0)
        return std::nullopt;

    return (v - n * vn) * eta - n * std::sqrt(delta);
}

/**
 * @brief 求解直线交点（如果不存在，则会返回最小二乘法解）
 *
 * @tparam N
 * @tparam T
 * @param[in] p1
 * @param[in] t1
 * @param[in] p2
 * @param[in] t2
 * @return std::optional<std::tuple<T, T>>
 */
template <std::size_t N, class T>
inline std::optional<std::tuple<T, T>> intersect_line_line(const vec<N, T> &p1, const vec<N, T> &t1,
                                                           const vec<N, T> &p2, const vec<N, T> &t2)
{
    mat<N, 2, T> A;
    mat<N, 1, T> b;
    for (std::size_t i = 0; i < N; i++)
    {
        A(i, 0) = t1[i];
        A(i, 1) = -t2[i];
        b[i] = p2[i] - p1[i];
    }

    auto ATA = A.transposed() * A;
    auto x = A.transposed() * b;
    bool invertible = math::solve(ATA.data(), x.data());
    if (invertible && (p1 + t1 * x[0]).distance(p2 + t2 * x[1]) < math::numeric::distance)
        return std::make_tuple(x[0], x[1]);

    return std::nullopt;
}

/**
 * @brief 获得平面点在三角形中的重心坐标
 * 分解 p = (1-u-v)A + uB + vC
 * uAB + vAC + PA = 0
 * [u v 1][ABx ACx PAx]^T = 0
 * [u v 1][ABy ACy PAy]^T = 0
 * u,v 是 右侧两个向量的叉积，并将 z 坐标归一化
 *
 * @tparam N
 * @tparam T
 * @param[in] p
 * @param[in] pts
 * @return std::optional<vec<3, T>> 重心坐标
 */
template <std::size_t N, class T> inline std::optional<vec<3, T>> barycentric(const vec<N, T> &p, const vec<N, T> *pts)
{
    vec<3, T> u1(pts[2][0] - pts[0][0], pts[1][0] - pts[0][0], pts[0][0] - p[0]);
    vec<3, T> u2(pts[2][1] - pts[0][1], pts[1][1] - pts[0][1], pts[0][1] - p[1]);
    auto u = cross(u1, u2);

    // 考虑退化情形
    if (std::abs(u.z()) < math::numeric::distance)
        return std::nullopt;
    return vec<3, T>{T{1} - (u.x() + u.y()) / u.z(), u.y() / u.z(), u.x() / u.z()};
}

/**
 * @brief 获得三角形上线性函数的梯度
 * 顶点函数值平移不影响梯度
 *      u1                u1 - u0
 *              =>
 *   u0    u2         0      u2 - u0
 *  grad 表示增加速度，沿着 v0v1 方向增加 u1 - u0
 *  (v0v1, grad) = u1 - u0
 *  (v0v2, grad) = u2 - u0
 *     (n, grad) = 0
 *
 *            [ v0v1 ]          [ u1 - u0 ]
 * A * grad = [ v0v2 ] * grad = [ u2 - u0 ]
 *            [  n   ]          [    0    ]
 *
 * @tparam N
 * @tparam T
 * @param pts
 * @param n
 * @return std::optional<mat<3, 3, T>>
 */
template <std::size_t N, class T> inline std::optional<mat<3, 3, T>> darboux(const vec<N, T> *pts, const vec<3, T> &n)
{
    mat<3, 3, T> A;
    A(0, 0) = pts[1][0] - pts[0][0];
    A(0, 1) = pts[1][1] - pts[0][1];
    A(0, 2) = pts[1][2] - pts[0][2];

    A(1, 0) = pts[2][0] - pts[0][0];
    A(1, 1) = pts[2][1] - pts[0][1];
    A(1, 2) = pts[2][2] - pts[0][2];

    A(2, 0) = n[0];
    A(2, 1) = n[1];
    A(2, 2) = n[2];

    return invert(A);
}

} // namespace geom

} // namespace xy
